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ABSTRACT 

Context. The revolution of helio- and asteroseismology provides access to the detailed properties of stellar interiors by studying the 
star’s oscillation modes. Among them, gravity (g) modes are formed by constructive interferences between progressive internal gravity 
waves (IGWs), propagating in stellar radiative zones. Our new 3D nonlinear simulations of the interior of a solar-like star allows us 
to study the excitation, propagation, and dissipation of these waves. 

Aims. The aim of this article is to clarify our understanding of the behavior of IGWs in a 3D radiative zone and to provide a clear 
overview of their properties. 

Methods. We use a method of frequency filtering that reveals the path of individual gravity waves of different frequencies in the 
radiative zone. 

Results. We are able to identify the region of propagation of different waves in 2D and 3D, to compare them to the linear raytracing 
theory and to distinguish between propagative and standing waves (g modes). We also show that the energy carried by waves is 
distributed in different planes in the sphere, depending on their azimuthal wave number. 

Conclusions. We are able to isolate individual IGWs from a complex spectrum and to study their propagation in space and time. In 
particular, we highlight in this paper the necessity of studying the propagation of waves in 3D spherical geometry, since the distribution 
of their energy is not equipartitioned in the sphere. 
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1. Introduction 

Internal gravity waves (IGWs) propagate in stably stratified 
fluids with gravity as a restoring force. They can be observed in 
laboratory experiments (e.g., Mowbray & Rarity 1967; Gostiaux 
et al. 2007; Perrard et al. 2013), planetary atmospheres, and 
oceans (e.g., Staquet & Sommeria 2002; Fritts & Alexander 
2003), and in our case, they are believed to exist in the radiative 
regions of stars. In solar-like stars, these waves are stochastically 
excited by turbulent convective motions in the outer layers, 
which leads to the formation of a rich spectrum. During their 
propagation, IGWs are known to transport angular momentum 
by radiative damping (e.g., Schatzman 1993; Zahn et al. 1997), 
corotation resonances (e.g., Booker & Bretherton 1967; Alvan 
et al. 2013), or nonlinear wave breaking (e.g., Barker & Ogilvie 
2010). They also affect the mixing of chemical elements in 
stars’ interiors (e.g., Press 1981; Garcia Lopez & Spruit 1991; 
Montalban, J. 1994; Charbonnel & Talon 2005). As a conse¬ 
quence, they are expected to influence the evolution of stars 
and particularly the evolution of their internal rotation profiles 
(Talon & Charbonnel 2005; Charbonnel et al. 2013; Mathis 
et al. 2013). In the solar case, IGWs are serious candidates 
for explaining the solid-body rotation of its radiative interior 
down to 0.2 R e (Kumar et al. 1999; Charbonnel & Talon 2005). 
1 They are also invoked to explain the transport of angular 
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1 Another hypothesis for explaining the solid body rotation of the 
solar radiative zone relies on the existence of a buried large-scale fos- 


momentum in sub- and redgiant stars (Talon & Charbonnel 
2008; Fuller et al. 2014) and in intermediate mass and massive 
stars (Pantillon et al. 2007; Lee et al. 2014). 

When progressive waves interfere constructively, they form 
standing modes called gravity (g) modes, whose measurement 
can provide information about the structure of the stars’ deep 
interiors. These individual g modes remain barely detectable 
at the surface of the Sun. Unlike acoustic (p) modes, mainly 
located near the surface, g modes are trapped in the inner 
radiative region, so we can accumulate a lot of information 
about the structure of this zone. Nevertheless, they are evanes¬ 
cent in convective regions and thus reach the surface with 
small amplitudes (Appourchaux et al. 2010, and references 
therein). As of today, the only reported detection for the Sun 
has been done in Garcia et al. (2007, 2008). They observed the 
asymptotic signature of g modes at the surface of the Sun, but 
we are not yet able to detect individual peaks (Garcia 2009). To 
improve our ability to detect them, we thus need to characterize 
their properties, time variabilities, and ability to tunnel through 
the solar convection zone better. 

Realistic numerical simulations help for understanding the 
properties of IGWs in stars. Their excitation by the pummeling 
of convective plumes at the top of the radiative zone has already 

sil magnetic field, whose existence and effect are still the subject of 
debates in the community (Gough & McIntyre 1998; Strugarek et al. 
2011; Acevedo-Arreguin et al. 2013). 
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been studied in 2D (Massaguer et al. 1984; Hurlburt et al. 
1986, 1994; Kiraga et al. 2003; Dintrans et al. 2005; Rogers & 
Glatzmaier 2005; Rogers et al. 2006) and in 3D (Saikia et al. 
2000; Brummell et al. 2002; Browning et al. 2004; Kiraga 
et al. 2005; Brun et al. 2011). Recently, Alvan et al. (2014) 
have shown the substantial progress made in high-performance 
computing, allowing a 3D solar-like star to be modeled in full 
spherical geometry from r = 0 to r = 0.97R©. They took the 
nonlinear coupling between the convective envelope, the radia¬ 
tive interior, and the dissipation along the wave propagation into 
account. Thanks to the use of a realistic density stratification 
in the radiative region, the properties and frequencies of the 
waves excited by the convection are those expected from theory. 
We seek here to characterize the spectrum excited better by 
doing a detailed analysis of its substructures. More precisely, 
we want to understand why we can see modes in the spectral 
space and not in the direct observation of the simulation (real 
space). Our guiding question concerns the distinction between 
progressive and standing waves and the transition between 
these two families. Finally, we see how the waves’ energy is 
distributed in the spherical star and what can be deduced from 
their excitation. 

In the present work, we use a method of frequency filter¬ 
ing that reveals the signature of IGWs in real space and that 
clarifies their behavior at different frequencies. In section 2, 
we show and discuss the complex and rich spectrum of IGWs 
excited by penetrative convection in our 3D model. In section 3, 
we describe our method of frequency filtering before presenting 
the associated results. We thus show the relation between our 
full 3D nonlinear simulations and the asymptotic linear raytrac¬ 
ing theory, often used to illustrate the propagation of internal 
waves in stars. We also show that waves of low frequency are 
attenuated while higher-frequency waves propagate up to the 
center and form modes. We finish by showing a property of 
internal waves that single-handedly warrants their study in 3D: 
at low rotation rate, the energy carried by IGWs is concentrated 
in planes of the sphere, whose inclination depends on the wave 
numbers. This property is broken at high rotation rates, leading 
to very complex 3D waves’ paths of propagation. 


2. 3D nonlinear simulations of IGWs 

We used the anelastic spherical harmonic (ASH) code (Clune 
et al. 1999; Brun et al. 2004) to solve the full set of 3D anelastic 
hydrodynamic equations in a solar-like star rotating at the ve¬ 
locity Oo = Do e z (Do = 414 nHz). The domain of computation 
includes both the radiative and the convective regions from r = 0 
up to r = 0.97 Rq. We provide a description of the ASH code and 
a presentation of the model used here in Appendix A. More de¬ 
tails can be found in Alvan et al. (2014). In this section, we focus 
on the study of the gravity wave spectrum excited in the radia¬ 
tive zone. This spectrum is rich and complex, and it varies as a 
function of depth. 


2.1. Wave pattern in physical space 

The use of a seismically calibrated ID solar model (Brun et al. 
2002) to initialize the 3D simulation ensures a realistic stratifi¬ 
cation in the radiative zone that allows us to study the proper¬ 
ties of the IGWs excited by the overshoot process quantitatively. 
To quantify this stratification, we define the Brunt-Vaisala fre¬ 


quency by 


N = 


IJJ 

\yP dr p dr J 


(1) 


where r is the radius, p, P, and g are the reference density, 
pressure, and acceleration of gravity, and y is the first adiabatic 
exponent (see Appendix A). In the whole paper, frequencies are 
plotted in Hz and not in rad/s. We represent the profile of N in 
Fig. 1. 



r/R 
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Fig. 1 . Profile of the Brunt-Vaisala frequency as a function of the 
normalized radius. Colored horizontal lines highlight the regions 
of propagation of the waves revealed by our filtering technique 
in Fig. 5. 


The linear asymptotic theory (e.g., Unno et al. 1989; Aerts et al. 
2010; Christensen-Dalsgaard 2011) predicts that the dispersion 
relation of IGWs is 


where co is the frequency of the wave (also expressed in 
Hz), and k the norm of the wavevector k = k r e r + kh whose 
radial and horizontal components are denoted k r and kh. This 
equation is satisfied when we neglect the effect of the Coriolis 
acceleration, which is possible here since the model rotates at 
the solar rotation rate, verifying 2D 0 = 2 x 414 nHz <$: co for the 
frequencies co of interest in this work. According to this relation, 
IGWs can propagate in the regions where N is real only, i.e. in 
radiative zones, and are evanescent in convective zones. Their 
maximum frequency is also limited by the maximum value 
max(V) = 0.466 mHz available for the current Sun (cf. Fig. 1). 

In Figs. 2(a) and 2(b), we show the radial velocity in a 
3D view and an equatorial slice of the simulated star with down¬ 
flows and upflows. We clearly see the convective envelope from 
about 0.71R© to the surface of the computational domain and the 
inner radiative zone. Since it is visible in Fig. 2(c), the velocity 
amplitude of the convective motions is much larger than the one 
of the waves (factor 10 3 to 10 5 ). Consequently, to visualize both 
convection pattern and IGWs in the top panels, we divided the 
radial velocity by its root mean square (rms) value at each radius. 


In the radiative region, we see what looks like concentric 
circles, which are in fact almost circular spirals. They corre¬ 
spond to the superposition of several wavefronts of gravity 
waves evolving at different frequencies. The comparison 
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Fig. 2. (a): 3D rendering of the simulated star. We represent the radial velocity divided by its rms value at each radius, (b): Equatorial 
slice. Like in the 3D view, IGWs pattern are visible in the inner radiative zone as quasi-circular spirals, (c): rms profiles of the total 
(solid line) and radial (dotted line) velocities as a function of the normalized radius, averaged over longitude, latitude, and time 
(about 10 convective overturning times), (d): Energy spectrum of gravity waves computed at ro = 0.5 R Q as a function of degree £ 
and frequency co. Ridges are formed by g modes of same radial order n, and we see that they tend to the maximum Brunt-Vaisala 
frequency at high order £ (dotted white line). The black solid line denotes the separation between g modes (above) and progressive 
waves (below). 
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Fig. 3. Percentage of the total energy (squared radial velocity) 
distributed in standing and progressive gravity waves. 

between these pattern and those predicted by the raytracing 
theory (see Sect. 3.2) shows that the IGWs visualized here 
have frequencies of about [0.01-0.05] mHz. This is quite low 
in comparison with the expected solar g modes’ frequencies 
(Christensen-Dalsgaard 2011). Deeper in the radiative region, 
the quasi-circular pattern changes into a more complex shape. 
We explain these observations in Sect. 3.3. 

2.2. Spectrum 

To describe the IGWs we can see into the radiative zone more 
precisely, we represent their spectrum so as to calculate their fre¬ 
quencies. It is obtained by applying a spherical harmonic trans¬ 
form, followed by a temporal Fourier transform, to the radial 
velocity v r (ro, 6, cp, t). For the moment, the radius ro is fixed, and 
we obtain the quantity v r (r 0 , l, m, oj), where t is the degree, m az¬ 
imuthal number, and oj the frequency of the wave. The quantity 
represented in Fig. 2(d) is obtained by adding all contributions 
in m quadratically, such as 

E(r 0 , (, oj) = ^ |v r (r 0 , £, m, a>)\ 2 . (3) 

m 

The length of the temporal sequence used here is about 30 days, 
which provides a frequency resolution of tu m i n « 0.38 x 10 -3 
mHz. The temporal sampling rate At = 1000s allows us to 
reach the maximum (Nyquist) frequency = 0-5 mHz and 
corresponds to ten times the time step of the simulation. 

The figure we obtained (Fig. 2(d)) is close to the spectrum pre¬ 
dicted by the linear theory (e.g., Provost & Berthomieu 1986; 
Christensen-Dalsgaard 2003). In the higher part of the spectrum 
(above the black curve), we see individual modes forming en¬ 
ergy peaks. Their frequencies correspond to g-mode frequencies 
calculated by the oscillation code ADIPLS 2 with very good ac¬ 
curacy (Alvan et al. 2014). Modes with the same radial order 
n - corresponding to the number of zeros of the radial eigen¬ 
functions - form ridges, particularly visible at high frequency. 
Moreover, as predicted by the dispersion relation (Eq. 2), the fre¬ 
quency upper limit corresponds to max(A0 = 0.466 mHz. The 

2 http://users-phys. au.dk/jcd/adipack.n 


lower limit of this “g modes” region is the effect of waves at¬ 
tenuation. It is shown here and can be understood as a cut-off 
frequency for the formation of g modes, under which the waves 
are sufficiently damped so that standing modes cannot form. We 
show in Appendix B that this curve behaves as [€(€ + 1)] 3/8 . 
Below this frontier, the spectrum is composed of progressive 
waves. When we compare spectra at different radii ro (Alvan 
et al. 2014), we see that the “progressive-wave region” decreases 
with increasing depth and finally disappears. We compare the 
proportion of energy distributed in both regions in Fig. 3. At 
r = 0.6R o , the whole energy is contained in progressive waves. 
Then, when we move deeper into the radiative zone, g modes ap¬ 
pear and progressive waves are damped. Below 0.3R o , the num¬ 
ber of progressive waves becomes negligible, and the whole en¬ 
ergy of the spectrum corresponds to standing modes. 

The richness of this spectrum shows that a wide frequency range 
is excited by convection. As a result, one may wonder why we 
only see one main pattern in real space (Fig. 2(a) and (b)). The 
raytracing theory predicts that IGWs propagate along different 
paths depending on their frequencies. Is there a way to visual¬ 
ize these paths in our simulations? By isolating some waves and 
visualizing them in the model, we show that we can learn more 
about the regions of propagation of IGWs in a full 3D geometry. 

3. Going further: frequency filtering of radiative 
zone 

Here, we have filtered our data by selecting only a narrow band 
of frequencies. In the following sections, we show that this 
method reveals important properties of IGWs that were not ac¬ 
cessible before. 

3.1. Method 

The idea of the frequency filtering can be discussed as follows. 
We illustrate the method in Fig. 4. 

Panel (A) shows a region of the simulated star belonging 
to the equatorial plane. We represent the radial velocity as a 
function of the radius and the longitude, right below the excita¬ 
tion zone (interface between convective and radiative regions). 
Several different gravity waves are propagating through this 
region, but they are placed on top of each other so we cannot see 
their paths clearly. Our aim is to separate the waves of different 
frequencies. To illustrate the manipulation, we have chosen two 
points (B and C) labeled by white crosses. Image (A) evolves 
with time and so does the signal at points B and C. 

Panels (B) and (C) represent the Fourier spectrum of these 
two points. For Point B, we see that the main peak is located in 
the blue shaded region, around 0.01-0.02 mHz. For Point C, the 
main peak is instead around 0.03 mHz (green shaded region) but 
we also see a secondary peak around 0.005 mHz (orange shaded 
region). This means that each point oscillates at the frequency 
of one or two waves passing through. 

Thus, to isolate a wave oscillating at frequency 0.02 mHz 
(for example), we calculated the Fourier spectrum of each point 
(i.e. N r x N e x Np = 569 x 512 x 1024 « 3 x 10 8 points), we 
applied a passband filter to the spectrum (a Gaussian centered at 
0.02 mHz), and we came back to the real space using an inverse 
Fourier transform. The width of the Gaussian filter is 2 x 10 -4 
mHz for this example. 
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0.02 mHz 0.03 mHz 0.04 mHz 



Longitude (rad) Longitude (rad) Longitude (rad) 


Fig. 4. Filtering of the same image at three different frequencies. Downflows are represented in blue, upflows in red tones. (A): 
selected zone of the star. (B) and (C) : frequency spectra of points B and C where we see one or two main peaks. The vertical axis 
represents the normalized radial velocity. The vertical red line indicates the position of the Gaussian filter (width = 2xlO~ 4 mHz). 
(Dl), (D2) and (D3): result of the filtering of image (A) at three frequencies. Wave beams are visible with an inclination that varies 
with the frequency, forming St Andrew’s crosses. The angles’ values are indicated in Table 1. 


The result is shown in Panels (Dl), (D2), and (D3) of Fig. 
4 for three different frequencies. For 0.02 mHz, we see that the 
isolated wave passes through Point B, and it is also the case for 
the wave at 0.03 mHz passing through Point C. We note a the 
angle between the direction of the wave beam and the gravity. 
This angle decreases from (Dl) to (D3), since we increase the 
central frequency of the filter. 


Filtered frequency 

N\ (mHz) 

<*th (deg) 

«m (deg) 

0.02 mHz 

0.03 mHz 

0.04 mHz 

0.15 + 0.02 
0.17 + 0.02 
0.18 + 0.02 

82.34 + 1.03 
79.84 + 1.21 
77.16+1.45 

85.4 + 2.29 
84.29 + 3.59 
82.2 + 1.79 


Table 1. Comparison between the St Andrew’s cross angles pre¬ 
dicted by the dispersion relation and measured in Fig. 4. Owing 
to the zooming effect, the a angles actually look smaller in the 
lower panels of Figure 4. 


This manipulation allows us to highlight one of the main prop¬ 
erties of IGWs. From the linear theory of ray tracing (Lighthill 
1978; Gough 1993; Christensen-Dalsgaard 2011), we know that 


IGWs propagate along beams whose paths are ruled by the dis¬ 
persion relation (Eq. 2). The shape formed by the beams at the 
point of the wave’s excitation is called a St Andrew’s cross 
(Lighthill 1986; Voisin 1991). The wave’s energy is radiated 
around an angle a± to the radial direction such that 

a th = arccos|^-j, (4) 

where N\ is the value of N in the region considered. This angle 
is visible in Panels (Dl), (D2), and (D3) of Fig. 4. To measure it, 
we calculate its tangent in each panel using the relation 

a m = arctan^), (5) 

where L is the length considered along the longitudinal direction 
and h along the radial direction. For the measurement of L, 
we have translated the longitude from radians to solar radii. 
The values found are indicated in Table 1 and compared to 
pseudo-theroretical values calculated from Eq. (4), where the 
value of N\ must be deduced from the figure. 

For a m , the main source of error is the measurement of 
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3D view 


Equatorial plane 





Fig. 5. 3D and 2D views of portions of the star filtered at different frequencies. On the right, ASH results are represented in gray 
tones, and we have superimposed the path of the corresponding wave calculated by the method of raytracing. Only the radiative 
zone is represented. For oj = 0.2 mHz, the green dotted circle represents the outer turning point, and the green arrows highlight 
the nonpropagation region. For oj = 0.3 mHz, the two blue dotted circles represent the outer and inner turning points, and the blue 
arrow shows the nonpropagation region. For the four frequencies, the propagation regions correspond to the colored horizontal lines 
in Fig. 1. 
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the distances h and L. The error shown here corresponds to the 
width of the beam. For a±, we also considered an error coming 
from the measure of N\ , which depends on the estimation of 
the starting point of the cross. Another possible source of error 
could be the quality of the filter, but since it has a very narrow 
bandwidth, this is negligible compared to the other factors. 

This comparison shows that we find a bias between both val¬ 
ues, where cc m is systematically larger. A future work will be 
dedicated to exploring more frequencies, in order to determine 
the extent of this difference precisely and to decide on the pre¬ 
cision of the linear dispersion relation used here. We have al¬ 
ready checked that the influence of the coriolis acceleration is 
too weak at this rotation rate to explain the observed difference, 
even though it goes in the right direction. 

3.2. Visualization of the high frequency wave pattern 

We now apply this method of filtering to a larger portion of 
the star to see if we can visualize other patterns than the quasi¬ 
circular spiral visible in Fig. 2 (a) and (b). Our results are shown 
in Fig. 5, in 3D and 2D views. To understand the results, we com¬ 
puted the paths in 2D of some IGWs using the method of ray¬ 
tracing. It consists in calculating the value of the wave’s phase 
along a path x(t) by resolving the Hamiltonian system 

( dxi _ dW 
d t dki ’ 

( 6 ) 

dki _ _dW 
< dt dxi ’ 

where W(x, k , t) = oj, and (x t , k t ) the cartesian coordinates of the 
position vector x and the wavevector k (e.g., Lignieres 2011). 
In our case, we employed polar coordinates, transforming the 
system (10) into 

dr _ dW _ 
dt dk r Vgr ’ 

dO _ 1 dW _ v g o 
dt r dk e r 

(7) 

dk r _ dW v g 0 
dt dr r 

dk e _ 1 dW Vgr 

- d 7 ~ ~~r~de~~ e ’ 

where v g = v gr e r + v g <9 ee is the group velocity of the wave at 
which the ray propagates. In this work, we use the raytracing 
theory in 2D, as employed by many other authors (e.g., Lignieres 
2011; Christensen-Dalsgaard 2011; Kosovichev 201 1), as a tool 
for understanding our simulations. When computing the ray 
paths shown in Fig. 5, we used the Brunt-Vaisala frequency 
profile presented in Fig. 1. Further developments of this method 
and its generalization in 3D will be the object of a future work 
(Mathis et al. 2015, in prep). 

The results shown in Fig. 5 illustrate the diversity of the 
paths followed by the different waves. The lower the frequency, 
the more the ray paths look like spirals. Another way to 
understand it is to say that we retrieve the cross shape at the 
point of departure of the waves, whose angle a with the radial 
direction decreases with increasing frequency. In particular, the 


yellow case (oj = 0.05 mHz) is the closest to the general pattern 
(Fig. 2 (a) and (b)). We thus confirm that what we see without 
filtering in the external region of the radiative zone is a sample 
of low-frequency waves. The other paths were not visible in Fig. 
2 (a) and (b) because the low-frequency part of the spectrum is 
dominant in amplitude (red tones in Fig. 2 (d)). By filtering out 
this part, we can visualize the region of propagation of the other 
waves, and especially of the resonant g modes. 

The dispersion relation (Eq. 2) imposes that the wave fre¬ 
quencies co are less than the Brunt-Vaisala frequency N. 
Consequently, we can define the limit of the propagation region 
by two turning points r in and r out such that N(r[ n ) = N(r out ) = oj. 
For instance, the wave oscillating at 0.3 mHz (bottom panel 
of Fig. 5) is confined more deeply in the radiative zone than 
the one at 0.05 mHz (top panel). This is also visible in Fig. 
1, where we indicate the regions of propagation of the waves 
shown in this section by colored lines. We note that the outer 
turning point is different from the penetration depth r ov , where 
convective plumes deposit the energy transferred to waves. We 
clarify the different radii defined here in Fig. 6. 



r/R @ 

Fig. 6. Representation of the overshoot region defined between 
r c , the first zero of the enthalpy flux, and the point r ov below 
which the amplitude of the flux drops to 10% of its most negative 
value. The radius rb cz corresponds to the point where the mean 
entropy gradient changes sign. The green vertical line represents 
the positions of r out for oj = 0.2 mHz. 


It is interesting to see that, although the excitation region is 
located at the base of the convection zone (penetration region), 
the perturbation is able to propagate in a “non-propagation 
region” (blue arrow and limits materialized by green and blue 
dotted circles). This phenomenon can explain the fact that 
high-frequency modes ([0.3-0.4] mHz) are less powerful in 
the spectrum (Fig. 2 (d)). In fact, a part of the initial energy 
transfered from the convective region to the waves is lost in the 
evanescent region, between the overshoot region and r out . 

A moving version of Fig. 5 shows that the paths visible 
here do not propagate. They oscillate at the chosen frequency 
(wavefronts propagate toward the surface at the phase velocity), 
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A co = 0.016 mHz 



Longitude (rad) Longitude (rad) 


Fig. 7. Distinction between progressive and standing waves. The top left panel is a zoom of Fig. 2(d) (same colorbar). Three 
horizontal dotted lines represent the frequencies chosen in Panels A, B, and C. In these panels, we show the radial velocity in the 
same region filtered at three different frequencies. Low-frequency waves are rapidly damped and cannot form g modes (Panel A). 
In contrast, for co = 0.1 mHz and co = 0.2 mHz (Panels B and C), we see modes oscillating between the two turning points (white 
lines). 


but the envelope remains stable. Thus, what we see here are 
stationary modes instead of progressive waves. In the following 
section, we show that the method of frequency filtering allows 
us to distinguish between these two families. 

3.3. Distinguishing between progressive and standing waves 

Along their propagation, IGWs are affected by a damping pro¬ 
cess that is proportional to the radiative diffusivity of the fluid. In 
the linear and asymptotic (co <$c N) approximations (Zahn et al. 
1997), the amplitude of a gravity wave propagating in a nonadi- 
abatic medium is damped by a factor e ~ T/2 where 

3 r° ut N 3 dr' 

T(r,£,o>) = [£(£+ 1)] 5 — (8) 

Jr u 4 r' 5 

We see that this depends on the radiative diffusivity coefficient k, 
but also on both the frequency co and the degree l. We introduce 
r out , the maximum radius, close to the external turning point 
r out (with r out < r out ), until when the JWKB approximation 
used to derive this expression can be assumed (see Fig. 6 and 


the detailed discussion in appendix B and in Alvan et al. (2013)). 

We know that g modes are formed by positive interferences 
between two progressive IGWs. This implies that the amplitude 
of these IGWs is high enough to reach the inner turning point 
n n and then to go back toward the surface (see Eq. 24). 

In Fig. 7 we show the result of different filtering of the 
same region of the star. This is a slice of the equatorial plane 
(6 = 7r/2) with the normalized radius as vertical axis and the 
longitude cp as horizontal axis. In the top righthand panel, we 
zoom in to the lower region of the spectrum presented in Fig. 2. 
We see that, at co - 0.016 mHz, the energy is fairly concentrated 
below the black line. The result is that the waves propagate 
over a finite distance before being completely damped. We draw 
the reader’s attention to the vertical axis that stops at 0.45 R Q 
since no wave is visible beyond. In contrast, for co = 0.1 mHz 
and co = 0.2 mHz, the main energy of the spectrum is above 
the black line (red tones), and in the real space, we clearly 
see that waves propagate between both turning points (white 
lines). They are excited at the top of the radiative zone, they 
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go from r out to where they bounce, and they come back to r out . 

We also note that the path of propagation changes between 
panels. In particular, the inclination of the rays steepens with 
higher frequency following the tendency discussed in Sect. 3.1. 
Moreover, we see that the number of rays increases with oj. We 
can define an effective wavenumber 4ff by 

4ff = (9) 

A 

where r ou t is the outer turning point, and A a wavelength that can 
be measured in the figures. We find 4ff ~ H for oj = O.lmHz 
and 4 ff « 14 for oj = 0.2mHz. 

The conclusion of this study is to confirm the intuition 
that the spectrum realized in our 3D simulations is made up of 
both standing modes and progressive waves, as one expects. 
Since the radiative damping is more efficient for high values 
of 4 the corresponding waves do not have enough energy to 
bounce at their inner turning point and to form modes. They 
thus stay simple progressive IGWs. Although not discussed here 
in detail, it is clear that viscosity will also damp the waves (in 
addition to radiative diffusion), since we have a Prandtl number 
close to one in the radiative zone (Vadas & Fritts 2005). 



Fig. 8. Schematic representation of the propagation planes of az¬ 
imuthal components of mode £ = 9. 


3.4. Energy concentrated in planes 


We now focus on a propagation property of IGWs in 3D. The 
results of this part highlight the importance of studying gravity 
waves in fully spherical geometry. The asymptotic theory of stel¬ 
lar oscillations (e.g., Gough 1993; Christensen-Dalsgaard 2003; 
Aerts et al. 2010) predicts that in the case where we can ne¬ 
glect the rotation of the star with respect to the waves frequen¬ 
cies (oj <$c 2Qo), the horizontal components of the wavevector 
verify 


( r sin Oky = m, 

k l = k l +k l = 


M + 1 ) 


r 2 


( 10 ) 


N 

X 


E 


o 

c 

0 

3 

cr 

0 



Using this property, Gough (1993) has shown theoretically that 0 5 10 15 20 25 

the ray paths of waves defined by (£, m) are confined in planes « 

forming an angle T 


sin 


€+ 1/2 


sin 


m 


y/€(€ + 1 ) 


( 11 ) 


with the polar plane. 3 These planes do not depend on the fre¬ 
quency of the waves, and they are independent of the dispersion 
relation. For instance, this means that this property is also true 
for acoustic waves. We show an example in Fig. 8, where we 
have represented the planes where modes (€ = 9, m = 0,3,5,9) 
are expected to be concentrated. We see that the highest values 
of m correspond to planes close to the equator. By definition, the 
m = 0 plane always contains the poles (axisymmetric case). 

We are going to test this theoretical result with our 3D nonlin¬ 
ear simulations. To illustrate it, we use the result of a filtering at 
frequency 0.3 mHz. The corresponding region in the spectrum is 
represented in Fig. 9. We show a zoom of Fig. 2(d) and indicate 

3 The exact expression found by Gough (1993) is the one using 
£+\/2, because he employed a general formalism that is different from 
the usual projection on spherical harmonics followed by the asymptotic 
approximation. Nevertheless, with respect to the precision of our re¬ 
sults, the approximation y/£(£ + 1) « £ + 1/2 is largely acceptable for 
the degree of the mode chosen here, £ - 9. 


Fig. 9. Zoom of the spectrum represented in Fig. 2(d). The dot¬ 
ted lines represent the bandwidth of the filter applied (Gaussian 
function centered at 0.3 mHz with full width at half maximum 
cr = 2 x 10 -3 mHz) that allows us to select the mode € = 9 in 
particular. 


the bandwidth of the filter. The mode located in the middle of 
this bandwith is £ = 9. We here recall that the spectrum shown 
includes all values of m (see Eq. 3). As a consequence, the red 
peaks of energy visible in this figure are formed by the super¬ 
position of all m components, whose frequency oj n ^ m is slightly 
shifted by the effect of the rotation. This rotational splitting is 
given by 


~ <*>nA 0 + 


m 


£(£ + 1 ) 


Go, 


( 12 ) 


where oj n z o is the central peak, corresponding to m = 0 (see 
Sect. 4.3.2 of Alvan et al. (2014)). 


Knowing the rotation rate of our model (1Q©) and the fre¬ 
quency considered here, the distance oj n f^ m - oj n j$ between 
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two peaks of different m and same £ is much less than the one 
between two modes of consecutive n values (see Alvan et al. 
(2014)). That is why we do not see any overlap in the peak and 
why the filter selects the whole azimuthal components of the 
mode. 

After having filtered the whole radiative zone, we need to access 
each plane of the sphere and to measure the corresponding 
dominant azimuthal number m. For example, to reach the plane 
09,5 = 31 (represented in green in Fig. 8 ), we apply a rotation 
matrix of angle 90 - 6 9,5 to the sphere, so that the desired plane 
is now located at the equator. We then apply a Fourier transform 
with respect to the longitude (p to obtain the spectrum in m. 
The result for four planes is given in Fig. 10. The white zones 
correspond to the dominant values of m. 

We see that these values are indeed different for each plane 
and that they correspond to those expected by the theory. In 
particular, the m = 9 peak is clearly dominant in the 69 $ plane 
and so is m = 0 in the 69 $ plane. That secondary peaks are 
also visible can be explained by the following arguments. First, 
for the top lefthand panel ( 69 $), we notice two weak peaks at 
m - 12 and m = 14. We assigned them to the mode £ = 14, 
which is partially taken into account by the filter (see Fig. 
9). Nevertheless, it is clear that the m = 9 mode is dominant. 
Then, the theory predicted that the existence of these planes 
assumes some approximation that is not verified by our model. 
The adiabatic approximation made by Gough (1993) is not 
true here since our waves are submitted to thermal and viscous 
diffusivities. This could result in a leakage of the energy from 
one plane to another. Moreover, the presence of rotation in our 
model is to be considered as a small disturbance. Indeed, by 
using our 3D ray tracing code (Mathis et al. 2015, in prep), we 
have shown that the propagation of IGWs in planes is no longer 
verified in the case of rapidly rotating stars (except for m = 0 ). 
Finally, we can invoke numerical noise because our spherical 
harmonic transform uses a finite number of mesh points, which 
results in a leakage of the energy from one pair {£, m) to another. 
Despite these limitations, this is a clear example of the need to 
model those waves in full spherical geometry. 

In Alvan et al. (2014), we showed that the energy of a 
given g mode tends to concentrate in high m components. We 
clarify this result in Fig. 11 by showing the distribution of energy 
as a function of latitude for the current model ( turb 2 ) and for a 
less turbulent model discussed in detail in that paper. Indeed, 
we see that for the model ref, the energy is concentrated around 
the equator. We suspect that this concentration of energy may 
be due to the shape of convection (banana cells) in low-latitude 
regions. If we increase the turbulence of the convection (model 
turb 2 ), the distribution becomes more homogeneous in latitude, 
but we still see a peak around the equator (colatitude 6 = 7r/2). 
For instance, this could imply a more efficient transport of 
angular momentum induced by the presence of IGWs in the 
equatorial regions. 


4. Conclusion 

The aim of this paper has been to understand gravity waves in 
solar-like stars. To this end, we analyzed our 3D simulations 
in order to deepen our understanding of the complexity of the 
spectrum excited. Here, we have explored and described the sub¬ 
structures of this spectrum by showing the co-existence of both 
progressive IGWs and standing modes. For the first time in sim¬ 


ulation, we were able to visualize individual IGWs in the inner 
regions of the stars. 

• Thanks to our frequency filtering method, we isolated some 
well-chosen waves in spectral space and shown that their en¬ 
ergy paths in the real space correspond to the one predicted 
by the linear raytracing theory. 

• We showed that it is possible to distinguish between g modes 
and progressive IGWs by measuring their damping rate. If 
the damping rate is too high, IGWs do not reach their in¬ 
ner turning point (whose depth depends on the frequency 
of the wave), and they cannot form resonant mode. In con¬ 
trast, we see in the figures that isolated waves taken in the 
high part of the spectrum (high frequency and/or high de¬ 
gree £) can bounce at r = and form a g mode. The cut-off 
frequency separating g modes and progressive IGWs scales 
with \£ {£ + 1)] 3/8 . This analysis gives us precise knowledge 
of the composition of the spectrum excited in our model. 
Consequently, it gives a rather good understanding of the 
same processes acting on the real wave spectrum of the Sun. 

• Our third result consists in the study of the spatial distribu¬ 
tion of the waves energy in the 3D radiative zone. We have 
shown that, according to the theoretical predictions of Gough 
(1993), waves of different degrees £ and azimuthal order m 
are distributed differently in space. Their energy is confined 
in planes, whose inclination with the polar direction varies 
with (£, m). The plane that is the closest to the equator is the 
one with m = £. Moreover, as noted in Alvan et al. (2014), 
the waves associated to high values of m have the highest 
amplitude, which may be explained by studying the reparti¬ 
tion of convective plumes. This result for the concentration 
of energy in planes is very important because it could guide 
our research into g modes at the surface of the Sun and our 
understanding of the angular momentum transport by gravity 
waves. 

It is important to highlight that, to perform these simulations, 
we adopted values of diffusivities that are necessarily much 
higher than their microscopic values. For this reason, some of 
the predictions of the model (e.g., concerning the energetical 
aspects of the spectrum of IGWs excited by convection) must 
be considered with caution. 

The direct perspectives of this work are to improve our analy¬ 
sis by taking the effect of rotation that modifies the ray paths 
and break the propagation in planes into account. This project 
requires improving both the simulation (exploring other param¬ 
eter regimes) and the theory. Indeed, our analysis is guided by 
the results of the raytracing theory but also by theoretical results 
about the impact of rotation on the excitation and propagation of 
IGWs (Belkacem et al. 2009; Mathis et al. 2014). We also plan 
to study the nonlinear aspects of IGWs, including for example 
triadic interactions or wave breaking (e.g., Staquet & Sommeria 
2002; Barker 2011; Bourget et al. 2013). Finally, the extension 
of this study to other types of stars is in progress and will provide 
a new source of asteroseismic predictions. 
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Fig. 11. Distribution of the energy transported by the IGWs as a function of the latitude for two different models. The vertical axis 
is normalized by the maximum value. 
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model, and v and k are effective diffusivities modeling momen¬ 
tum and heat transport by subgrid-scale (SGS) motions that are 
unresolved by the simulation. Their profiles are represented in 
Fig. 12 . 


Appendix A: ASH code and 3D model description 


The ASH code solves the set of 3D anelastic hydrodynamic 
equations in a spherical geometry. The equations are fully non¬ 
linear in velocity, and we linearize the thermodynamic variables 
with respect to a spherically symmetric and evolving mean state. 
We denote p, P, T , and S as the reference density, pressure, tem¬ 
perature, and specific entropy. Fluctuations in this reference state 
are denoted by p, P, T, and S . We assume a linearized equation 
of state 


p _ P T P S 

p P T yP c p ' 


(13) 


and the zeroth-order ideal gas law 


P = RpT , 


(14) 



r/R 

o 


where y is the first adiabatic exponent, c p the specific heat per 
unit mass at constant pressure, and R the gas constant. 

We also introduce the local velocity v = (v r ,V 0 ,Vy>) ex¬ 

pressed in spherical coordinates (r, 0, y) in the frame rotating at 
constant angular velocity Oo = Do e z . The set of hydrodynamic 
equations solved by ASH in the present case is 

' (a) V. (pv) = 0, 

(b) P^ +(v.V)vj = -pVS-p^g-2pQ 0 xv 

- V.£>- [VP-pg], 

dS 

(c) pT— +pfv.v(s + 5) = pe + V. [K r pc p V (t + ?) 

+KpTVS + K ( ,pf VS j + 2 pv [eijeij - 1/3 (V.v) 2 ], 

(15) 

where equation (a) is the continuity equation in the anelastic 
approximation, (b) the momentum equation, and (c) the energy 
equation. 

We define the gravitational acceleration g , the viscous stress ten¬ 
sor D defined by 

Dij = -2 pv(e u - 1/3 (V.v)tf y ), (16) 

with eij = 1/2 (djVi + diVj) the strain rate tensor and <5 (/ the 
Kronecker symbol. The bracketed term on the righthand side of 
the momentum equation is initially zero since the system begins 
in hydrostatic balance. For this equation, we used the Lantz- 
Braginsky-Roberts (LBR) form (e.g., Lantz 1992; Braginsky 
& Roberts 1995) advocated by Brown et al. (2012), who have 
shown that the classical anelastic formulation may introduce a 
bias in the conservation of energy, in the case of strongly sta¬ 
bly stratified atmospheres. In this formulation, we use the re¬ 
duced pressure To - P/p instead of the fluctuating pressure P 
and neglect the buoyancy term pcoV (s/G 3 ) associated with the 
non-adiabatic reference state in the radiative zone. Thus, the new 
buoyancy term pgS/c p only contains the contribution of entropy 
fluctuations, and the contribution due to pressure perturbations is 
included in the reduced pressure gradient. In the energy equation 
(c), K r is the radiative diffusivity based on a ID solar structure 


Fig. 12. Radial profiles of the effective diffusivities k (thermal) 
and v (viscous). The model used here is called turb2 in Alvan 
etal. (2014). 

The diffusivity kq is also part of the SGS treatment in the 
convective zone and is chosen such that the entropy flux carries 
the solar flux outward at the top of the domain. It is negligible 
in the radiative zone (Miesch et al. 2000). Finally, the volume 
heating term pe represents the energy generation by nuclear 
burning with e = £o T k . The constant £o is set such that the 
radially integrated heating term equals the solar luminosity at 
the base of the convection zone, and k-9 allows us to take both 
contributions of the p-p chains and CNO cycles into account 
(Brun et al. 2002). 

The computational domain extends from r^t = 0 to r top = 
0.97 R q where R Q = 6.9599 x 10 10 cm is the solar radius. The 
boundary conditions at the top of the domain are torque-free 
velocity conditions and constant heat flux (Brun et al. 2011). 
The inner boundary conditions are chosen to deal with the 
central singularity at rb 0 t = 0. Conditions on the poloidal and 
toroidal components of the mass flux impose that only € = 1 
mode can go through the center, and thermal condition are also 
adapted. This treatment is explained in detail in Alvan et al. 
(2014). 

In the model presented here, the spatial resolution is 
Af r x Nq x N(p = 1581 x512x 1024, with a non-uniform 
radial grid chosen to resolve the turbulent convection, the 
oscillations of IGWs in the radiative zone, and the steep drop 
of diffusivities at the interface between both zones (Alvan et al. 
2014). 

The 3D simulation is initialized using a reference state de¬ 
rived from a ID solar structure model (Brun et al. 2002). The 
set of equations (15) is then solved in expanding the velocity 
and thermodynamic variables in spherical harmonics Y^ m (G, p>) 
for their horizontal structure and using a fourth-order finite- 
difference approach on a nonuniform grid for their radial part. 
For the time integration, we use an explicit Adams-Bashforth 
scheme for the advection and Coriolis terms and a semi-implicit 
Crank-Nicolson treatment for the diffusive and buoyancy terms 
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(Glatzmaier 1984; Clune et al. 1999). The reference state is 
updated during the simulation with the spherically averaged 
perturbation fields. 

After a few convective overturning times, a balance is es¬ 
tablished between the contributions to the total flux of the 
different physical processes (Brun et al. 2011; Alvan et al. 
2014). It allows us to see that the convective flux becomes neg¬ 
ative at the base of the convective zone, which is the signature 
of the overshoot process (e.g., Zahn 1991), where downflows 
penetrate the radiative zone slightly owing to their inertia. This 
induces a perturbation in velocity and temperature at the top of 
the radiative zone, which propagates in the form of a gravity 
wave. With the bulk excitation process (Lecoanet & Quataert 
2013), overshoot is responsible for the excitation of IGWs in 
solar-like stars. 


Appendix B: Understanding the cut-off frequency 
between modes and progressive waves 

The goal of this appendix is to derive the dependence of the cut¬ 
off frequency that separates modes and progressive waves, co c , 
on their latitudinal degree {€). 

We develop the radial component of the velocity on spherical 
harmonics: 


u r 


Z u n(,m ( r ) Y hn (i9, <p) exp (iojt). 


(17) 


Next, we get from the system formed by the linearized momen¬ 
tum, continuity, and heat transport equations in the adiabatic 
case 

d 2y ¥ tm (r) 7 

+ k 2 v . { (r) (r) = 0, (18) 

_1 ? 

where = p 2 r z u rAm , and 


k 


2 _ 

V\l ~ 


IN 2 \£(£ + 1) 

U 2 ) r 2 


(19) 


From then on, we focus on the low-frequency asymptotic regime 
in which co <$c N. We follow Zahn et al. (1997) by applying the 
JWKB and the quasi-adiabatic approximations that allows us to 
derive the expression for u r .^ m : 


u r ,(,m = Ae^r ip P| m| (cos0) 


X cos cot + 


f ?OUt \ 

kvjdr exp 


r (r, G (o) 


, ( 20 ) 


where A^ m is an amplitude coefficient, which includes the nor¬ 
malization of spherical harmonics, and P™ are the associated 
Legendre polynomials. We recall the radiative damping expres¬ 
sion (Eq. 8) 


3 r° ut N 3 dr' 

r(rJ,co) = [W+l)P k ——. ( 21 ) 

Jr OJ 4 r' 5 

Since the JWKB approximation falls at the two turning points 
(rin, r out ), we introduce the two radius (n n , r out ) (with r m < r m 
and r out < r out ) between which it can be applied (see the detailed 
discussion in Berry & Mount (1972) and Alvan et al. (2013)). 


Then, we derive the horizontally averaged kinetic energy of 
IGWs using results derived by Zahn et al. (1997): 

(r) = \p { U lm) g <p = («L m) B „ > (22) 

where ( mmm )o t(p = l/4n sinOdOdcp. Using Eq. (20), it becomes 

En-,e,m (r ) = ([4 m| ( cos 6 tf) g r ” 3 ex P [ _T kr, t, oj)] . 

(23) 

We then define the cut-off frequency that separates modes and 
progressive waves with the following criterion: 

EK;£,m (hn) = KE K .^ m (r out ) . (24) 

When K <$c 1, the wave has lost a large part of its kinetic energy 
and cannot reflect to form a standing mode. Using Eq. (23), it 
becomes 


r r N ( r °ut) r out 

exp [-T (r in , i, w)] = K — — - = a. (25) 

N(r m )r-i 

Assuming that r m and f out vary weakly with co c for co c N, we 
finally obtain 


f 


rout dr' 

kN — 

r' 3 3 

[W+IW 


In o' 


= p[({l+ 1)]» (26) 


as observed in power spectrum of our direct 3D nonlinear ASH 
simulations (Fig. 2). 
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